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Based on a continuum meclianical model for single-layer graphene we propose and analyze a 
microscopic mechanism for dissipation in nanoelectromechanical graphene resonators. We find that 
coupling between flexural modes and in-plane phonons leads to linear and nonlinear damping of out- 
of-plane vibrations. By tuning external parameters such as bias and ac voltages, one can cross over 
from a linear to a nonlinear-damping dominated regime. We discuss the behavior of the effective 
quality factor in this context. 



I. INTRODUCTION 



Advances in fabrication and detection techniques 
have enabled a wide range of experimental realiza- 
tions of carbon-based nanoelectromechanical (NEM) 
resonators^ However, to optimize their operation, 
an increased understanding of dissipation mechanisms 
is needed. For NEM resonators in general, several 
processes leading to linear damping (LD) have been 
investigatecPl'^. Specifically for graphene, at high tem- 
peratures, ohmic losses in the metallic gate and the 
graphene sheet have been argued to limit the quality 
factoi'^. Recently, the focus has sh ifted to study quantum 
aspects of mechanical motiorP^liiJ, such as mechanical cat 
state^I^!^ which require a more detailed understanding of 
dissipation and decoherence mechanisms. 

Since graphene-based resonators exhibit nonlinear be- 
havior, one can expect the damping also to be amplitude 
dependenl!i2HlIl, Nonlinear damping (NLD) was reported 
in recent experiments on graphene and carbon nanotube 
resonators^. However, little is known about the underly- 
ing physical mechanism, and typically phenomenological 
models are employecP^^t^^. In these models, the resonator 
is coupled to a bath of harmonic oscillators. For couplings 
that depend quadratically o n the re sonator amplitude, it 
is known that NLD emerge^i^ilfiHII. 

For carbon-based resonators such a coupling naturally 
arises if the strain couples linearly to the degrees of free- 
dom of some subsystem, which can be regarded as a bath. 
Two ex ample s are the interaction between phonons and 
electron j^^l^^l and the coupling of mechanical modes. The 
relative importance of the two mechanisms is a priori not 
known and will also depend on the details of the experi- 
mental realization. 

In order to quantify the importance of the mechani- 
cal dissipation channel for NLD, we analyze the coupling 
between flexural modes and in-plane phonons. We show 
that it leads to a quadratic coupling and, consequently, 
to both LD and NLD. Whether LD or NLD dominates 
is determined by the ratio of vibrational amplitude and 
static deflection. We give an estimate for the expected 
crossover between LD and NLD, which can be experi- 
mentally verified. 




FIG. 1. (Color online) Schematic view of a suspended 
graphene membrane over a trench in an insulating substrate. 
A metallic gate is used for actuating the resonator. In-plane 
phonons are created in the suspended region and dissipate 
energy as they propagate away. 



II. MODEL AND METHOD 

We consider a graphene sheet of length L and breadth 
6, suspended over a trench of width £ (cf. Fig. [l]). The 
van der Waals attraction between the graphene and the 
substrate clamps down the sheet outside the suspended 
regiorp21122] xhe trench is modeled by allowing the sheet 
to freely displace vertically in this region. Since out-of- 
plane displacement is accompanied by in-plane stretch- 
ing or compression, flexural motion is converted into in- 
plane phonons in the suspended region. The clamping 
constrains the out-of-plane motion over the substrate, 
but still allows for small in-plane displacements. Con- 
sequently, in-plane phonons created in the suspended re- 
gion transport energy away from this region. In contrast 
to a phenomenological modeling approach we can relate 
dissipation to specific properties of the substrate and the 
graphene-substrate coupling. These properties can be 
obtained independently by theoretical or experimental 
means. 

The dynamics of graphene NEM-resonators are well 
described by the continuum theory of 2D-membrane^2ll, 
For a resonator made from a sheet lying in the a;?/-plane, 
this theory is conveniently formulated in terms of the in- 
plane displacement fields u(x,y),v{x,y) in the x— and 
y— directions, respectively, and the displacement field 
in the 2;— direction, w{x,y). The equations of motion 
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follow from the free energy = J dxdy + J-"s] where 
J^b — f lAwp is the free energy density associated with 
pure bending and J^s = j '^ij^ij associated with 
stretching of the membrane. The symmetric 2D strain 
and stress tensors are here defined as 



The free energy ^ leads to a coupling between flexural 
vibrations and in-plane motion via the coupling energy 
^coup = (7i/2)u which is nonlinear in the flexural 

vibration amplitu de. T his coupling leads to NLD of the 
flexural vibrationJSElHIZI. 



wl/2, (la) 



and 



(^xx — (Ag + '2^J-g)^xx + ^G^yy, <^ xy = '^fJ-G^xy, 

^yy — ('^G + 2^G)eaa + ^G^xx , 



(lb) 



respectively. Spatial derivatives are denoted by sub- 
scripts, i.e., = du/dx. The expression for the free 
energy^ which is similar to that for large deflections of 
a plate^, contains three material parameters, the bend- 
ing energy k « 1.1 — 1.6 eV, and the Lame parameters, 
/iG « 146 N/m and Ag « 48 N/m for graphenepHHI. To 
study qualitatively the effect of phonon radiation into the 
supporting substrate, we assume for simplicity a quasi ID 
situation where variations in y— direction are disregarded. 
This would be valid for a wide sheet where deviations 
from this assumption is confined to the regions around 
the edges. In this case we have only the displacement 
fields u{x,t) and w{x,i). In any realistic functioning de- 
vice, there is some small amount of built in strain. In 
practice, this implies that the energy contribution from 
the bending energy is always negligible for the lowest ly- 
ing flexural mode^23 Hence, to a good approximation 
we have for the quasi ID graphene resonator attached to 
a substrate the free energy density 
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+ ]^K{x){u~usf +E,^t[w], (2) 

where we have defined Ti = Aq + 2/iG- The potential 
foxti^'j] accounts for interactions used to actuate the res- 
onator. The second to last term couples the graphene 
displacement to the substrate displacement us{x,y) in 
a harmonic approximatiorP^, which largely allows us to 
obtain an analytical description. 

The function K{x) restricts this coupling to the sup- 
ported region, i.e., K{x) = KoQ{\x\ — £/2) with Q being 
the Heaviside step function. The substrate is modeled 
as an elastic half-space and displacement at the surface, 
s{x,z = 0,t) = {us,v s,ws), is given in terms of a re- 
sponse functioiP^EHEI]^ 

Sp(f, z = 0,0;) = y" ^-^i?^^(f - x',a;) 

X a„z{x',u}) . (3) 

Consistent with the ID model of the graphene sheet, only 
us{x) = ^^^i^dy us{x,y) is considered. Within the har- 
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monic approximation, 



A. Equations of motion 

The equations of motion for the out-of-plane and in- 
plane vibrations resulting from Eq. ([2| are 



T d 

PGW ~ Ydx (^"■^'"'^ + ^^^) ^ -^dc + /ac COs{nt) 



T d 

Pqu - ^^"'"^ ^ ^^^^ ^ ~K{x) (u - us/b) , 



(4a) 

), 
(4b) 



where fdc{x) and fs^dx) cos{D,t) are the static and time 
dependent parts of the actuation force. Typically, elec- 
trostatic actuation is used, resulting from a time de- 
pendent back-gate voltage of the form Vhg{t) = Vdc + 
V'accos(rit) with Vdc ^ V"ac- To simplify the analysis, we 
assume the equilibrium stress field resulting from /^c to 
be spatially uniform and equal to the tensile stress Tq 
on the boundarji^H. Generally, at a given back-gate bias 
voltage, the resonance frequency f2o(Vdc) depends on ini- 
tial stress and contains a shift due to electrostatic forces. 
This so-called tuning behavior will be further discussed 

in Sec. mra 



Since Eq. (4b) is linear in u, the influence of the en- 



vironment can be accounted for by a Green's function 
embedding technique. The solution. 



u(x, t) 



dx 



dt'G{x, X , t - Y ^■'",a;'^(a;', t') 



K{x) {u - us). 



. (5) 

is given in terms of the in-plane response function G, 
which contains information about the attachment to the 
substrate via Eq. Th e speed of sound in graphene is 
denoted by c = ^rTi / pQ , where pQ is the mass density 
of graphene. 



B. Flexural mode dynamics 

Next, we consider the fundamental flexural mode and 
set w{x,t) = q{t)(f){x) for \x\ < £/2 and zero otherwise. 
The mode shape (p is normalized to the length of the res- 
onator. Upon projecting Eq. ( [4a] ) onto the fundamental 
mode, an ordinary differential equation for the vibration 
amplitude q is obtained. Further, moving to a rotating 
frame, we write q{t) = [qo + ^ (gi(t)e*"* + ql{t)e''"*')] 
and q{t) = [gi(i)e'^" - g5;(t)e"*^"] . Inserting these 
expressions into the equation of motion and performing 
the averaging yields an equation for the slowly varying 
amplitude qi |13j . which contains memory terms related 
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to linear and non-linear damping. As the time-scales for 
flexural motion and in-plane phonons are well separated 
(f2o ^ cji), the memory terms can be eliminated. This 
procedure corresponds to a Markov approximation^. It 
is convenient to define new quantities 
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dx 
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dx 
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where G{x,x',u}) = (27r)-i/ dTG(a;, x', r)e*"^ is the 
Fourier transform of the in-plane response function. 

We obtain an equation of motion for the complex en- 
velope function 



by the ratio of the overlap integrals defined in Eq. (|6|, 
which are purely geometrical quantities, and the ratio 
between the vibrational amplitude and the static deflec- 
tion. For a small static deflection, it is therefore expected 
that NLD dominates the damping caused by phonon ra- 
diation. Similarly, the dimensionless ratio 



V = 



V ^0 



(10) 



measures the relative importance of the two nonlineari- 
ties in Eq. ([7p^. For fj < the well-known bifurcation 
of the Duffing equation is present, while for fj > y/3 this 
bifurcation vanishes. The ratio fj is also a purely geo- 
metrical factor, apart from the weak dependence of fig 
on the static deformation of the graphene. 
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For finite temperatures this equation has to be sup- 
plemented by noise forces, satisfying the fluctuation- 
dissipation relations. The thermally induced vibrations 
can l ead t o an additional broadening of the response 
curve^i^'^. In order to obtain a lower bound of LD and 
NLD we will work in the limit of zero temperature. In 
Eq. Q, the coefficients m ~ Pg^^, a, 7 and r] denote 
the suspended mass, the Duffing elastic constant, linear 
and non-linear damping, respectively. They are given in 
terms of x ^ follows 
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x(o) + ^xm 
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q'o^lmxin) , 
2Im xi^fl) . 



(8a) 
(8b) 
(8c) 



Here, the bare Duffing constant is given by ao — 
(T16/2) J dx (t),x(x)'^- The driving strength is g — 
^ dx (j){x)f^c{x). In accordance with our previous sim- 
plifications, we neglect the small polaronic shift of fio, 
which is proportional to Re Xi s-iid additional shift of 
a due to the broken symmetry in the presence of static 
deflection. Equation ([t]) is similar to the equa tions used 
to model NLD in micromechanical resonatorJSEH and 
recent experiments on carbon-based resonatorJ^, the dif- 
ference being the dependence of the damping coefficients 
in Eq. ([s]) on the driving frequency. 

In Eq. ([7| the prevailing damping mechanism is deter- 
mined by the ratio 



5 = 



mi\ 



Im x{2^l) \qT 



47 8Im x{^) 



Here, jg™'^'^! denotes the maximum amplitude of the re- 
sponse for a given driving strength. Thus, 5 is determined 



Numerical method 



To compute the overlap integrals (|6| we first consider 
the Fourier transformed response of the substrate ([3| 



us(x,a;) 



L/2 fc/2 fc/2 

dy' I dy 
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{x — x' ,ll>) axz{x' ,uj) . 
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In the second step, in order to get a purely ID response 
function, we have approximated the y'-dependence of 
Cxzix' , y') by the mean value laxz and defined Rxx{x — 



x',^)=\ /!{/2 dy' J^'' dy R^^ix 



6/2 
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'X',y-y',uj^. The 
response function i?^^ for an elastic half-space is known 
analyticalljEilSSIlI] a,n(j niainly depends on the longitudi- 
nal and transversal sound velocities of the substrate (see 
Appendix |A| . 

Evaluating Eq. (11) at discrete positions {xi}^ leads 



to the linear system 

Kug(a;) = - [I - KM(a;)]" 



^(cj)Ku(w) , (12) 



which can be solved for ua^{xi, uj). Here bold-face symbols 
denote vectors of length TV, e.g., u = [u{xi), . . . , u{xm)\ 
and double struck symbols are N x N matrices. In 



particular, lij 
(27r)-2i?^^(a;, - 



= S. 



Ky = K{xi)5i_j and — 
Using this result and the dis- 



cretized version of the equation of motion ( 4b I one ob- 



tains an equation for the in-plane response function 



(9) (j^^. — G{xi,Xj,uj) 



-ujH - c^L + 
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PG 



[I - KM(cj)] 



(13) 
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where L is the discrete second derivativd^. Approximat- 
ing the integrations in Eq. Q by numerical quadratures, 
one finally obtains 



graphene and substrate parameters 



xm = -*'G(-r!)* 



(14) 



with $i = 'Ix'^'^xlj.-x ' which allows the computation of 
X for a given geometry. The parameters entering the 
equation of motion can then be calculated using Eqs. 
([S]). Following Ref. HH we set 7 — ^/{mCla), 7) = -qQa/a, 

9 = ffv^^/^O' ^ = ^I^Q, and q ^/ajm^. In the 
limit of weak LD, 7^1, the response of the resonator 
is determined solely by the dimensionless parameters 77, 
g and f2, describing the nonlinear damping, the driving 
strength and the driving frequency. 



III. RESULTS 



graphene mass density 
Ac, + 2/iG 
Si02 mass density 
Si02 sound velocities 

coupling strength 



PG 7.6 X 10"^ kgm^2 

Ti 340 N m"^ 

PS 2.2 X 10^ kg m"^ 

cl/c 0.28 

ct/c 0.18 

Ko 1.82 • 10^°N m"^ 



resonator parameters 



total length 

length 

width 

distance to gate 
tensile stress 



L 
I 
b 
d 

To 



2 p,m 
1 jj,m 
1 fim 
330 nm 
0.34 N m-i 



TABLE L Graphene and resonator parameters used for the 
calculations in Figs. |3] and [4] Graphene and substrate param- 
eters are taken from Refs . I35l and 1361 



To quantify the influence of LD and NLD, we con- 
sider the setup shown in Fig. [T] with a back-gate volt- 
age Vbg = Vdc + 14cC0s(rit). The fundamental-mode 
shape is taken to be 0(a;) = cos(7ra;/^), which gives 
ao = iTxT^^h I {M^) . Within a parallel plate model for 
electrostatic actuation, the force acting on the graphene 
sheet is given by 



2{d + q{t)<i>{xyf 



{yl, + 2ydcKc cos(r!i)) , 

(15) 



where G(yj) = eo/((i + w) is the capacitance of a par- 
allel plate capacitor with plates being separated by the 
distance d -\- w and eg is the vacuum permittivity. The 
distance is determined by the depth d of the trench and 
the flexural displacement w of the resonator. In the sec- 
ond line we further assumed Vdc ^ V'ac , which is typically 
found in experiments. The force can be separated into a 
static and a time-dependent part, / = /dc + /ac cos(rit) 
with /dc oc and /ac oc VdcKc, respectively. Since 
the displacement, which is on the order of a few nanome- 
ters, is much smaller than the trench depth, the force can 
be expanded in powers of w. Accordingly, the driving 
strength in Eq. Q becomes g — 2\/2i?6eoKicKc/(7i'c^^)- 
Moreover, the static displacement can be found by solv- 
(4a) and ( |4b[ ) in the static limit (see Appendix 
_ This yields q^ w ypie €qVIJ {'k'^ d^T^) . Note the de- 
pendence on the tensile stress Tq; go becomes smaller for 
increasing tensile stress. 

In the following, we consider a graphene resonator 
with dimensions and parameters as given in Tab. |l] We 
checked that the results do not change, for larger values 
of the total length L. Using Eqs. Q and (10) we obtain 
a/ao ~ 0.64 and ij Ri 7 ■ 10^^. The latter implies bi- 
stable behavior of the resonator. In general, these values 
depend sensitively on the geometry of the graphene sheet 



ing Eqs 
B 



and on the substrate. Our results provide a "best case" 
estimate, since the substrate is treated as a semi-infinite 
medium and the trench is modeled by the position de- 
pendent coupling K{x). Lifting these restrictions will 
lead to a stronger response of the substrate, and more 
dissipation. 



A. Resonance frequency 



As described in Sec. |II A| the resonance frequency 
^^o(Kic) depends on the initial stress and the bias volt- 
age. The dependence of fio on bias voltage, the so called 
tuning curve, is a characteristic feature of NEMS devices. 
It is a result of the competition between softening (de- 
creasing ilo) due to the electrostatic force [Eq. (15)], and 
stiffening (increasing fio) due to the Duffing nonlinearity 
of the graphene sheet. 

To obtain the tuning curve, we separate static and dy- 
namic contributions to the displacement fields. 



w{x, t) 
u{x, t) 



W(3{x) 
Uq{x) - 



- 5w{x, t) , 
5u{x, t) 



(16a) 
(16b) 



and insert these expressions into the equations of motion 
given by Eqs. (|4|. The static solutions, wq and mq, are 
calculated in Appendix [B] Further, we expand the static 
force fdc(x) up to first order in 5w, 



■''''^ 2id + wo^ ^ {d + wo)3 



(17) 



The resonance frequency is then obtained by collecting 
terms, which are linear in the vibration amplitude Sw. 
There are three such terms, which contribute to the res- 
onance frequency. 



C!g(Vdc) = f^^(0)+Ar!Lch. 
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FIG. 2. (Color online) Resonance frequency Jlo vs. bias volt- 
age. Symbols denote results of numerical calculation. The 
dashed (red) and dashed-dotted (blue) lines show the contri- 
butions of mechanical stiffening and electrostatic softening for 



To = 10 Ti, respectively. Parameters are given in Tab.^ 
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Og(0) = 



Af^Lch.(14c) 



Al^e\(^dc) 



7o^ 
PgP ' 



3m 



(18b) 
(18c) 
(18d) 



The three contributions are due to initial strain, mechan- 
ical stiffening and electrostatic softening, respectively. 
Since the static deflection go depends on the bias volt- 
age Vdc, the last two terms yield the voltage dependent 
tuning behavior. 

Figure [2] shows the tuning curve for the parameters 
given in in Tab. |T] For voltages, Vdc > 10 V, the reso- 
nance frequency (squared) is mainly determined by the 
mechanical stiffening, which scales with 



V^* while the 



(18) 



softening term scales with V^^ according Eqs. 

Depending on the specific geometry and the initial 
stress, the resonance frequency of the resonator may be 
substantially tuned using the bias voltage. Since the lin- 
ear and nonlinear damping coefficients given by Eqs. ([8| 
depend on frequency, the magnitude of LD and NLD will, 
in principle, also be influenced by the tuning curve. In 
order to disentangle the influence of rio(^c) and the cou- 
pling to the in-plane phonons, we will only consider a con- 
stant resonance frequency fio = f^o(O) = •JTqJpq,{ti / 1) 
in the following discussions (see Appendix O for the in- 
fluence of the tuning on the quality factor). 
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FIG. 3. (Color online) Ratio 5 of nonlinear (NLD) and linear 
(LD) damping terms according to Eq. (|9|; a) bias voltage and 
b) ac voltage dependence. The thin dashed and dashed-dotted 
lines show the asymptotic behavior for strong LD and NLD. 
A crossover between the two regimes is achieved by changing 
the bias voltage. Parameters are given in Tab. |T] 



B. Damping ratio 

The relative importance of LD and NLD, which is 
quantifled by 5 deflned in Eq. ([9]), is determined by 
the ratios Im x(2^^)/(8Im x{Vl)) and Igf^'^'^l/go- The 
former weakly depends on the geometric details. For 
small f2 one can expand Im x(^) odd powers of O. 
As Im X is proportional to the density of states of the 
substrate phonons, Z3(i7) oc f2, we expect on symmetry- 
grounds for a quasi-lD geometry, that Im x(ri) oc Vi^ . 
Consistent with this expectation, we obtain numerically 
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Im x(2fl)/(8Ini x(n)) w 0.93. 

The maximum amplitude qf^^^ can easily be found 
from Eq. ([7| in the steady-state limit, which yields an im- 
plicit equation for the magnitude l^ij of the steady-state 
amplitudelA. Sweeping the driving frequency, the max- 
imum amplitude is attained when d\qi\/dft = 0, which 
results in the cubic equation 
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45=|g— |(47 + 77|gr'^P). 



(19) 



Here, 7 and g depend on the bias voltage Vdc via qo and 
/ac, respectively. However, note that only g depends on 
the ac voltage. Due to the different dependencies of go 
and \q™'^^ \ on the bias voltage, one can achieve a crossover 
from NLD to LD dominated behavior by increasing the 
bias voltage. This is shown in Fig. [3^. In the limit of 



smaU Vdc, |grl ~ {"ig/vy" oc Fa/V^/' and 6 > 1, i.e., 
NLD dominates. For large Vdc, Igf""''! ~ 5/7 « VacV^j"^ 
and S goes to zero with increasing Vdc- Since the static 
displacement is determined only by the geometry and the 
bias voltage, and the maximal amplitude additionally de- 
pends on the ac voltage, the crossover can also be realized 
by tuning T4c, which is shown in Fig. [3]d. Equating the 
expressions for in the two limits gives an estimate 

for the crossover for both voltages. Additionally, due to 
the dependencies of qo oc Tq^ and f^o oc \/Tq on the 
initial tension Tq one finds that the damping ratio S in- 
creases with increasing tension in both regimes {S oc Tq 
and J oc Tq in the LD and NLD regime, respectively). 
Thus, the non-linear damping is enhanced for larger Tq. 



C. Quality factor 

To quantify the energy loss we consider the quality 
factor Q = no{E±) / {E±) , which measures the time- 
averaged dissipated energy {Ej_) normalized to the av- 
erage energy {Ej_) in the flexural modes. The nonlinear- 
ities render Q amplitude dependent. To get a worst case 
estimate, we use the maximal amplitude. In the slow 
envelope approximation we find 



^0 + 
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(20) 



The nature of the damping influences Q. In the LD dom- 
inated regime, S <^ 1, Q is independent of the vibrational 
amplitude, Qld ~ 'm^o/l- In contrast, for 6 > 1 one gets 
Qnld ~ 4ml7o/(»7k™'"'P) for V > 1. Thus, Q increases 
with decreasing driving strength. This agrees with the 
conclusions of Ref. 01 

Figure shows the quality factor as a function of bias 
voltage for constant Vac- As expected, Q decreases with 
increasing bias and excitation voltages and its behav- 
ior with regard to applied voltage changes qualitatively 
at the crossover between LD and NLD regimes. The 
asymptotic LD behavior limits the maximally attainable 
Q-factor, which is indicated by the gray area. We also 
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FIG. 4. (Color online) Quality factor Q vs. bias voltage, a) 
Q calculated from Eq. ( |20[ | and b) with additional voltage 
independent damping, = +Qo^ with Qo = 10^. The 
gray area indicates the region of attainable Q-factors. The 
dashed lines correspond to the behavior in a). Parameters 
are given in Tab. |T] 



compare to the case where the LD is additionally caused 
by a mechanism that does not depend on the bias volt- 
age leading to Qp. In this case the effective Q-factor, 

Fig. 13 

factors. The qualitative difference between LD and NLD 
is still present and should be experimentally observable. 
Most importantly, by decreasing Vac the maximally at- 
tainable Q-factor, which is determined by other damping 
mechanisms can be approached. 



i ^ +Qo^, has a cutoff for small Vdc as shown in 
which further limits the region of attainable Q- 



IV. CONCLUSIONS 

In conclusion, we have studied coupling between flex- 
ural vibrations and in-plane displacements as a physical 
mechanism for damping of flexural modes in graphene 
resonators. A characteristic consequence, which influ- 
ences the behavior of the dependence of the quality fac- 
tor on bias and excitation voltages, is the competition 
between static deflection and vibrational amplitude. We 
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note that the same type of behavior would naturahy oc- 
cur for any dissipative process which couples linearly to 
the strain; for example, Ohmic dissipation induced by 
synthetic gauge fieldJ^^. The cross-over should allow 
for an experimental verification of this class of damping 
mechanisms. 
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Appendix A: Response of an elastic half-space 

The displacement response at the surface of an elastic 
half-space to a stress acting on the surface is given in 
terms of a response function by Eq. ^ . If the stress is 
directed parallel to the x-axis, the spatial Fourier trans- 
form of Eq. ([3]) reads 



(Al) 



where k = {kx,ky) is the surface wave vector. The re- 
sponse func tion R ^xik, uj) for finite frequencies is explic- 
itly given bjJsnill] 



Rxx{k,Lo) 
with 

PL,T(w,fc) : 

S{OJ,k) : 



kl 



PT(^,fc) kl : 

^ S{uJ,k) 4 fc2 ^ pT,{uj,k)k'^ j 

(A2a) 



UJ 
UJ 

cl,t 



f ie — fc2 
2fc2 



(A2b) 



Ak'^pi,{uj,k)pT{uj,k) , 
(A2c) 



where cl and ct are the longitudinal and transversal 
speeds of sound, respectively, and the infinitesimal e > 
ensures causality. Notice that pl,t and S{uj, k) depend 

only on the modulus k of the wave vector k. The response 
function in real space is then 



Rxx{x,uj)= I d^kRxx{k,uj)e"'-'' 
2iTi / d 



2 \ f,J-i^^v) + Qyhi^^y) 



(A3a) 



Here, we defined 



-Hy2 \ctJ J S{u},k) 



dk 



\/a;2 +y'^ J PT(w,fc) 



Ji{kyG!^Ty^ 



(A3b) 



(A3c) 



where Ji is a first order Bessel function of the first kind. 
Note, that 

Ix{x,-y) = Ix{x,y) , Iy{x,-y) = -Iy{x,y) . (A4) 



The expressions given in Eqs. (A3) are a very conve- 



nient starting point for the numerical evaluation of the 
response function used in Sec. |IIC[ 

The zero-frequency response can be directly calculated 
in real spac^^. One finds 

_ , _ 1 2(4 ~ cl)x^ - cjy^ 
4:Trpsc^ (cl - c^){x^ + y^y/^ 



Appendix B: Static displacement 

In the static limit, the equations for the in-plane and 
out-of-plane displacements Q within the suspended re- 
gion become 



Tiu^xx + yC*. {w%) =0 , (Bla) 
dx [{2u^x + W%) W^x] ^fdc , (Bib) 



with vanishing boundary conditions a,t x = ±£/2 for the 
out-of-plane displacement. To find the proper boundary 
conditions for the in-plane displacement, we need to con- 
sider the coupling to the substrate in the non suspended 
region. Here, the equation for the in-plane displacement 
( 4b ) is given by 



Tiu. 



K{x){u{x) -us/b) ^ 



(B2) 



Following the same line of reasoning as in the main text, 
the static substrate response can be written as 



L/2 



us{x) 



-L/2 



dx' 



^Rxx{x-x')Q{\x'\-l/2)h{x') . 



(B3) 



with h{x) ~ KQiu—ug/h) and Rxx{x—x') being the static 
response function for an elastic half space given by Eq. 
( A5 ) integrated over y. To treat the problem analytically. 



we convert Eqs. (B2| and (B3) into a local equation for 



the in-plane displacement. In the limit of very strong 
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coupling to the substrate, the spatial variation of h{x) is 
small, in which case 



L/2 



dx' — 

R.Ux-x')Q{\x'\-l/2)h{x') 



(2^) 



-L/2 



h{x) 



L/2 



-L/2 



dx' 



■R^^{x-x')Q{\x'\-il2) . (B4) 



This makes it possible to solve for h{x) in terms of the 
in-plane displacement u{x). One finds 



h{x) = 



l-Ro(x)Ko 



u{x) , 



(B5) 



L/2 



where i?o(a;) = (27r)-2 / dx' R^^{x - x')Q{\x'\~ £12). 

-L/2 

This expression is valid outside the suspended region and 
is approximately given by h{x) « —1/Rq{x), which as- 
sumes KoRo{x) 3> 1. Consequently, the equation for the 
in-plane displacement, Eq. (B2|, is modified to become 



(B6) 



for |a;| > £/2. Thus, the effect of the substrate is re- 
duced to that of a spring with a spatially varying spring 
constant. The displacement u is expected to decay ex- 
ponentially to zero in the clamped region with a decay 

length A = ^ Rq{x)Ti. For the substrate parameters 
given in Table |lj this amounts to A w 100 nm. As a 
consequence, within a distance of 100 nm from the edge 
of the suspended region the in-plane displacement u{x) 
is essentially zero. To a good approximation, we there- 
fore assume vanishing boundary conditions for in-plane 
displacement at |a;| = £/2. 

Setting u(x) = (To/Ti)x+Au{x), where the first terms 
accounts for initial strain in the graphene, the bound- 
ary conditions are w{x = ±^/2) = and Au{x = 
±£/2) = 0. Using the Ansatz w{x) — qo(f>{x) with 
4i{x) = -\/2 cos TTx/^, the in-plane equation ( Bla| ) reads 



Alt . 



90^ 



(B7) 



Consequently, the in-plane displacement will be given by 



Au{x) 



£2 



dx' sin^ : 



2P 



QqX 



(B8) 



Inserting this expression into Eq. (Bib) and we obtain 



90 I iprTo 



Tiql 



2^/2 



dc 



(B9) 



This is a purely algebraic equation for the static deflec- 
tion. In the limit qo <^ ^\f^ ~ "^"^ ^^"^ ^ ~ ^1^'^ 



10 : 



o 



■a 



10 : 



td 10' 



10' 



E V 7^ = 510 

~- n T = 10'^ T 

r 1^ iu 

A r = 10"^ T ^' 



-I— r-pTTJT 
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FIG. 5. Static deflection go vs. bias voltage for thre e values of 
initial tension. The linear approximation [Eq. ( BIO l] is shown 



as dashed lines in the figure, while the squares and triangles 
correspond to the full numerical solution of the static problem. 



and To/Ti = 10 ^,the cubic term can be neglected and 
90 « /dc- 

To compute goj we need to consider the electrostatic 
interaction with the back gate. The static force acting 
on the graphene is given by Eq. (15). Considering the 



limit go ^ '^i we obtain for the static displacement 



9o 



\^Tod^ 



(BIO) 



III In Fig. 



which is the expression given in Sec. 
linear approximation (dashed line), given by Eq. ( |B10| 
is compared to the full numerical solution of Eq. ( |B1[ ) 
(squares and triangles), which takes the substrate into 
account. The linear approximation remains valid in the 
displayed interval for the two larger values of initial strain 
To /Ti , while a more significant deviation is apparent for 
the lowest value of the strain. 



Appendix C: Influence of tuning and initial tension 
on the quality factor 



In Sec. |III A| we discussed the voltage dependence of 
the resonance frequency (tuning curve) and showed that 
the frequency can be substantially tuned by changing the 
bias voltage Vdc- Since the linear and nonlinear damp- 
ing constants given by Eqs. ([s]) depend on frequency, the 
quality factor will also depend on the tuning. In order to 
quantify the influence of the voltage dependence of the 
resonance frequency on Q, Fig. [6] shows the quality factor 
for constant flo — rio(O) (dashed lines) and flo{Vdc) (full 
lines). One sees that deviations between these two cases 
appear only for larger voltages (Vdc > 20 V). Moreover, 
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FIG. 6. (Color online) Quality factor Q vs. bias voltage cal- 
culated from Eq. Q for Vac = 10"* V. The full and dashed 
lines show the result for a voltage dependent rio(Vdc) and con- 
stant Qq = f2o(0), respectively. Parameters are given in Tab. 

m 



the qualitative behavior and the cross-over from NLD 
to LD behavior remains unchanged. This confirms our 
statement in Sec. |III C[ that the behavior of Q is dom- 
inated by the damping coefficients 7 and 77 rather than 
the voltage dependence of fio- 



Additionally, Fig. |6] shows the quality factor for a 
smaller value of the initial tension. In this case, the qual- 
ity factor is decreased for all values of the static bias volt- 
age. In the limit of large LD, this is due the increased 
static deflection (see Eq. (BIO)). In the opposite limit, 
the quality factor is independent of the static deflection, 
and the decrease in quality factor is instead a result of the 
decreasing resonance frequency rio(O) oc y/To. Further- 
more, as argued at the end of Sec. |III B[ the cross-over 
between NLD and LD is shifted toward lower values of 
the bias voltage, signifying a decrease in the importance 
of NLD for lower tension. 
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